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A ROUGH EARTH SCATTERING MODEL FOR 
MULTIPATH PREDICTION 


1. INTRODUCTION. 


The most important phenomena to be considered in a model of radio wave communica- 
tion between earth satellites are scattering from the surface of the earth. In Section 2 
we describe in detail the model we have derived implemented on a computer to 
predict the field received after reflection from a rough, spherical earth. In Section 2. 1 
our discussion follows that of Beckman and Spizzichino^ but differs from theirs and 
from all others. Section 2. 2, in two important respects. First, the scattering integrals 
are not approximated, but are computed numerically; second, the domain of integration 
is the appropriate region on the surface of the earth. The work of Beckman «nri others 
deals with the flat earth approximation which is inaccurate for satellite to satellite 
communication. 

The results of the calculations are discussed in Section 3. Calculations have been per- 
formed at VHF frequencies and for terrain which could be described as marshy land. 
Various different altitudes are considered as well as several values of the physical 
parameters which describe the roughness of the earth. The results show tha t: 

a. Rough surface scattering calculations must be performed over 
a spherical earth when satellites are involved. 

b. There is a definite dependence on the values of the roughness, 
a, and the correlation length, T. The ratio of o/T can not be 
taken as an arbitrary small number as has been done previously 
in some calculations. ^ 


> 
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Continued. 


For some types of pulsed co mmu nication systems, the multipath 
effect is largest for incident angles (the angle of the incoming 
ray measured from the normal to the surface) which are approxi- 
mately 80 degrees. Some previous analysis have considered the 
worst case to be at an incident angle of zero degrees. 
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2 . THE SCATTERING OF RADIO WAVES FRO M THE SURFACE QZ 

THE EARTH. 

[ 1 ] . 

In this section we first summarize the work of Beckman and Spizzichino which 
describes a method for obtaining the average value and the mean square value of the 
amplitude of an electric field received at a point after reflection from a rough surface. 
It is assumed that the wave illuminating the surface is a plane wave and that the sur- 
face is a plane except for random irregularities which are normally distributed. This 
theory is then applied to the problem of a spherical wave impinging upon the surface 
of the earth by dividing the illuminated portion of the surface into small elements 

such that: 

a. Each element of surface 6S is large compared with the hori- 
zontal dimensions of the irregularities 

b. Each element of surface is sufficiently small so that the incident 
angle of the incoming wave does not vary appreciably from point 
to point on that element, i.e. , the spherical wave may be appraxi 
mated by a plane wave over each element. 

The contributions from each element of surface are added together to yield the total 
field reflected from the surface of the earth. 

Section 2. 1 summarizes Chapters 3 and 5 of Beckman's book; the application to scat- 
tering from thfl surface of the earth is made in Section 2.2. 


i 
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2.1 Scattering From Elemental Areas. 

Let a surface S be given by 


z = t(x.y) 

where x, y, and z are cartesian coordinates. We will consider £ a random process of 
the space variables x y with zero mean, and the statistics of the scattered field will 
be determined by taking certain averages over this random surface. 


In order to obtain expressions for the average amplitude ami the average mean square 
value of the field scattered up to the receiver, we require a knowledge of the first and 
second order statistics of the random surface, i.e. , the distribution function 

F i(0 = p («<*•» £ 0 

where P denotes the probability measure, and the second order joint distribution 

f 2 ( c r0 = p ( c (*r y i) * c i' c (V y 2 ) £ « 2 ) 

(The dependence of F g on the points ( x r y x ) ( x 2 » y 2 ) haS been ^PP*® 88 ®^) Th ® 
probability density functions are defined by 
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2.1 — Continued. 


and 



* r a M«) 
*2 


We assume that the surface may be adequately described by the normal distribution 
(although other distributions have been used); we obtain therefore 



where a 2 is the variance, and C = C(r) is the correlation coefficient, r being the 
distance between y^ and ^x^, y g)* 
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2.1 — Continued. 

We choose 


C(r) 


2 

T 



where T is called the correlation distance. 

Let E be the field incident on the surface £(x, y), and let its amplitude be a 
1 

plane wave: 




r 


where the time factor e 1 wt has been suppressed, 


k 

t 2 « 1 

1 = * k i 


is the propagation vector, X denotes wavelength, and 

r = (x,y,z) 

is the radius vector, as indicated in Figure 2-1. 
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2.1 — Continued. 

The angles ft, ft, and 0 are indicated in Figure 2-2. k is in the direction of 
1/2 3 ^ 2 

reflection at the surface, and = \k^ J = k^. 

Let P be the point of observation, and let R ' be the distance from P to a point 
x, y, C (x.y) 011 the surface S. Then the amplitude E of the scattered field is 
given, by 

J j(*%- *£) d ° 

S 


where 


* = 


ik 2 R' 

e 





and E and dE/dn are the field and its normal derivative on S. If P is at a great 
distance from S, then we can make the far field approximation 


k ft R ' » k R 
2 2 0 


k • 
2 


r 


where R Q is the distance from the origin to P (see Figure 2-1). The field on the 
surface and its normal derivative are found by means of the Kirchhoff approximation 
(B(ll), B(12), section 3.1). [B(ll) refers to Equation 11 of Reference 1.] 
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2 


Figure 2-2 „ Definition of Scattering Angles 
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2.1 — Continued. 


E = (1 +R) E 2 


and 

6E -* -* 

JE = i(1 ’ R ) E 1 k 1 # n 

where n is the normal to the surface z = £ (x, y), and R is the reflection coefficient 
of a smooth plane which depends on the incident angle, the reflecting material, and 
the polarization of E^. The reflection coefficients were originally derived by Fresnel 
and are often referred to as the Fresnel coefficients. See Section 2. 3 for a discussion 
of the Fresnel coefficients. 

Applying these approximations to Equation (3) we obtain 


E 2 ( p ) 


ikR 
i e 

0 


0 


//■ 


(Rv - p) • n 


S 



< 4 ) 


where 



p " k l + k 2 


1 
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2.1 — Continued. 


It is convenient to normalize E 2 by E 2Qt the amplitude of the field reflected by a smooth, 
perfectly conducting plane. Setting R = land0 1 = 0 2 in Equation (4) yields 


ikA cos 0 


20 


*R, 


1 ikR n 
- e 0 


where A is the area of the projection of S onto the x, y plane. We obtain 


P 



1 

A cos 0 



. -* -» 
iv • r 


e 


dx dy 


(cf. B(l), Section 3. 2), where 

v = k[sinQ l -sin0 2 cos0 3 , - sin 0 g sin 0g, - cos 0 X - cos 0 2 J 
a = (1 -R) sin 0 + (1 +R) sin 0 2 cos 0 g 
b = (1 4R) cos 0 2 - (1 -R) cos 6 1 
c = (1 +R) sin 0 sin 0 g 

ds = V l"+ J! 2 + C 2 dx dy 
* y 

and the region of integration has been approximated by a rectangle in the x, y plane. 
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2. 1 — Continued. 


In order to make progress with this formulation, it is necessary to treat a, b, and c 
as constants, in which case, (cf. B(10) Section 3.2 and B(7) Section 5.4) 


P 



Y X 



-* -♦ 
iv • r 
e 


dx dy 


(5) 


where 



1 + cos 0^ cos 0 2 - sin 6^ sin 0 2 cos © 3 
cos 0J ^cos Q 1 +cob 0^ ~ 


<R> - average of R over the surface 


<r> = ^P^C'JRtC'JdC' 

we approximate <R>«r( 0 = 0jJ. 

The average value of pis found to be 


<R> F 

<t» = —r 



ixv 

x 


♦ Syv 
y 


iCv 

<e > dx dy 


= *0 *( V z) 


( 6 ) 
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2.1 


— Continued. 


where 


*0 = 


<R> F sin v X sin v Y 
3 x y 


iCv. 


( T z) 9 < e Z> ‘ J 


iv £ 


e 


P x (0 dC 


is the characteristic function of the probability density P, . 

1 


For the normal distribution we have (B(3), Section 5. 3) 


x (\) * 


e 


, 2 2 
-1/2 o v 

z 


( 7 ) 


The average power received is proportioned to pp*; therefore we require an approxi- 
mation of i op* s 


<flP*> = 


<|rI 2 >f 3 2 y y x x 


L L L L 


[ iV x( X 1 - X 2) +1V y( y i' y 2)j 


• *2 ( V \) "i *2 dy ! “ y 2 


( 8 ) 
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2.1 


—Continued, 


where 


X 2 <u, v) = <e 1 


ivC, iuC, 


°? ® iv iuC 2 

e e 


•/ / 


p 2 («x*f 2 ) d fi d < 2 


For the probability density given by Equation (2), we have (B(40), Appendix C) 


X 2 ▼) 



♦ 2Cuv + v 


2 


)] 




so that 



e 



(1-C) v b 


1 


(») 


It is convenient because of convergence problems of certain integrals below to calculate 


D(p) e <pp*> - <p> <p*> 


2-1 2 



ESL-PR53 


2. 1 — Continued. 


D(p) 


<|r| >f. 


Y Y X X 

//// 

-Y -Y -X -X 


[ lv x ( X 1 ’*2) + lv y 


[x 2 (V - V 0 • X ( V *) X *( V z)] dx l d *2 dy l dy 2 


In order to simplify this integral, Beckman introduces two transformations, 
is defined by 


r i = X 1 “*2 


U 1 ~ X 1 + X 2 



U 2" y i 


♦ y. 


whereupon Equation (10) becomes 


D(p) 


<|r| 2 >f 3 2 




iv v. * iv r_ 
x 1 y 2 


[ X 2 (V -\) ' X ( V .) X * ( V z)] dr l dT 2 


X - y 2>] 

( 10 ) 

The first 
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2,1 — Continued. 

(Certain terms have been neglected here, Beckman justifies this by c l a iming that the 
quantity in brackets is very small except for and x^ near zero. ) The second trans' 
formation is: 

r 1 = r cos <p 
t 2 = rsin <p 

which yields 


D(p) 


<|r| 2 >f 3 2 

A 



iv r cos p + iv r sin <p 

ex y 


• [X 2 ( V Z- V 2 )‘ X ( V I ) X ( V 2 i )] Td, ’ dT 


where now the region extends over the entire xy plane, the justification being that the 
quantity in brackets is very small for large x. The formula 



ix cos <p + iy sin <p 


e 


dip 


then gives 


D(p) 
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2. 1 , — Continued. 

I 2 2 ~ 

where v = Wv + v 

xy yf x y 

(There is a misprint in B(33), Section 5.2; A should be A. Also <|r| > appears in our 
formula because of B(8), Section 5.4.) The function in brackets is, from Equation (7) and 
Equation (9) 

> 

Mv v *)' x ( v i) X *(M " 

and 

2 2 T2 no 

g = v z ° " hr 


e ' -e 


(cos + cos 0 2 ) j 


At this point previous authors (Beckman, Section 5 . 3) have taken an approximation to 
Equation (11) which is valid only when g is very large. Since all values of g may be 
impor tant in a given problem, we choose to evaluate Equation (1 1) numerically. 


2.2 Scattering From the Surface of the Earth. 

The results of the previous section will now be applied to obtain approximations of the 
field at the receiver after the reflection of a spherical wave from the surface of the earth. 
We the total surface of reflection to be broken up into surface elements 6 S, each 

of which is approximated by a section of a plane. Each of these plane surfaces reflects a 
certain amount of energy to the receiver according to the results of the previous section. 
These amounts are added to yield the total energy reflected into the receiver. 
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2.2 — Continued. 

The reflecting surface is indicated in Figure 2-3. 


Recall that 


D<p) = <pp*> - <p> <p*> 


and 



where 

ikR 0 

OA e cos 0 , 

E = — I 

20 nR Q 


and 0 A is the area of the surface element 6 S. 


Defining 


D ( E 2 ) - < E 2 V> ' < E 2> <E 2*> 


we have 

d ( e 2 ) - l E 20 l 2 d <p> 


> 
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2.2 — Continued. 


which from Equation (11) is 

D(E 2 ) = 0 1 ) <|H|> 2 

00 

. 2» f 3 2 J J 0 ( v xy r )[ X 2 ' XX *] TdT <12> 

0 



where we have multiplied by the divergence coefficient D^, r^, 6^ (cf. B(8X 
Section 11. 3) to take into account the curvature at the element 6 S (r and r are 

1 A 

the distances from the transmitter and receiver) and by the shadowing function, 

which will be discussed in Section 2.4. The shadowing function accounts for the 
shadow cast by one part of the surface onto another part of the surface. 

) * <V< E 2>* 

which, from Equation (12), is approximately 


The total power at the receiver is then <E E* 

• 2 2 


> = 


1 


D (E 
i V 2 


< E 2 E 2 > 


■ ? ffrfr ** 9 ») 


GO 

<|R | 2> / J o(V) [ x 2 ' XX *] 

0 


XX* rdrd S 


+ <e 2 > <e 2 > 


(13) 
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2.2 — Continued. 

This three-fold integral must be evaluated numerically. For this purpose it is con- 
venient to introduce the coordinate system indicated in Figure 2-4. B is the specular 
point (the point at which 6^ = 9^, 0 g = 0), and P is an arbitrary point in the scattering 

region, A and C are the earth subpoints for the transmitter T and the receiver R. We 

measure P by the angles 6 and <p which are obtained as follows. Construct the great 

circle through P which is perpendicular to the great circle ABC at M. Define 0 = BM 

and (p = MP. The region of integration is the intersection of the two spherical caps 

defined by those points P for which 0 < 0 , where cos 0 = RE/RE + H and y<y, 

o o o 

where cos y = RE/RE + HI. From the identities of spherical trigonometry, it is 
o 

easily shown that 


dS = cos d 0 d <p 

In order to find <E 2 E 2 *> (s® 6 E< l uat *on (13)) we add <E 2 > <E 2 *>to D^E,^. Since 

<E_>is nearly zero except for specular reflection, we may determine <E >by modifying 
2 2 

Ej by the path loss together with an appropriate reflection coefficient. 

This reflection coefficient is given by (cf. B(4) Section 12.3). 

R s = <*s> DR 

where R is the reflection coefficient of a smooth plane earth, (see Section 1. 3), D is 

the divergency coefficient which takes account of the earth's curvature, and p is the 

s 

value of at the specular point, which is given by exp -1 / 2 8, from Equation (6). 
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CENTER 
OF EARTH 


Figure 2-4. Coordinates System and Associated Angles and Distances 
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2.2 — Continued. 

We have included also a shadowing function (cf. Section 2.4) so that 

<E > <E *> = 8 ~ t , D 2 H 2 S(8,) 2 

2 2 (R + S) 2 ° V V 

where R and S are the distances of the transmitter and receiver from the specular 
point (see Figure 2-4). 

2. 3 Fresnel Reflection Coefficients. 


f21 

The Fresnel coefficients which are used in the calculation are the exact expressions. 
The reflection coefficient for the component of the wave with electric vector perpendic- 
ular to the plane of incidence is given by 



("2 q ~ "x a 2 C ° B 9 l) 8 * "2 p2 

(“ 2 t,+M l a 2 008 e i )* l ‘2 p2 


where 0^ is the incident angle from the normal, and the subscript 2 refers to air and 1 
refers to the reflecting surface material, p and q are found from the relations 


2 

q 




2 2 
a 2 8in °1 


2 

q 









1/2 
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2.3 — Continued. 

and a and /3 are given by 



with the 1 or 2 subscript specifying which value of p, c and a to use. p, c and a are 
the permittivity, permeability and conductivity in rationalized MKS units. The phase 
change 6 X is given by 


tan 


2 «1 M 2 a p cos e. 



For the E field component parallel to the plane of incidence 



and the phase change 6 n is given by 
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2.3 — Continued. 


In the case when a^/c « 1 the reflection coefficients reduce to 




( K el - K e2 8l ° 2 9 x) V2 - 


K _ cos 0., 
e2 1 


K - K sin 0 
V el e2 


\V2 

l) + V^2° OSe i 


K „ cos 
el 


Rii « 


'»■ V s * ( 


k l/2 


K - K „ sin 0, 
el e2 li 


K , cos 
el 


>1 * ^2 ( K , 


a/2 


. , - K _ sin 0, 
el e2 ly 


where 


K = — 
e e_ 


and c is the permeability of free space. This approximation is used in Beckman and 

Jl] 


Spizzachino 1 J with K e replaced by 


K = - s - + i 60 X a 
e <f 

o 

f21 

For the approximation with K = c/e . Stratton 1 shows that reasonable values are 

e o 

ob tain ed for R ± and R|| for frequencies above 1 MHz and below 10 4 Mhz. The exact 
expressions for R x and R t | are valid at any frequency. 
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2.4 Shadowing Function. 

The shadowing function, S (0^, is included in the calculation of <EE*> and <E><E>* 
to account for the screening of parts of the surface by other parts as shown in Figure 2-5. 
The shaded area in Figure 2-5 does not contribute to the scattering area and thus any 
function which contains a factor of area, A, should be modified to contain A. 

s(0i) is defined as 


S 0 1 = <S (x y)> 


where S (x y) is one at a point which is illuminated and zero at a point which is not 
illuminated. 

We use a shadowing function derived by Beckman. ^ For a normally distributed 
surface 8(0^ is given by 

- ^ tan 0^ erfc (k cot 0J u 

■ft)-* L J 


where 



T = correlation length 
a = RMS value of the roughness 
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2.5 Discussion of Approximations. 

Several approximations are present in the formulation we have presented. Here we 
discuss these and indicate areas where refinements might be made in the future . 

a. The Kirchhoff approximation was applied to the integral repre- 
sentation (Equation (3)) of the field E^P) to obtain the field E 
on the surface. This requires that the radius of curvature of 
the surface (which is approximately o/T) is small compared 

to a wavelength; hence, the approximation breaks down if the 
surface includes sharp edges. Brekhovskilch (B(2) Section 3.3) 
gives as a criterion 4nr cos 0 where r is the radius of 

curvature, and 0 the local angle of incidence. 

b. The point of observation must be far from the scattering surface. 
This approximation was invoked to obtain k R ' in the expression 
for ij) in Equation (3). For satellite communications this approxi- 
mation will always be satisfactory. 

c. The Fresnel coefficient and the function F g (0^ 0^, 6^ ) appearing 
in Equation (4) have been replaced by their average value and 
removed from under the integral sign in Equation (5). 

d. Two approximations have been introduced in the transformation of 
the integral Equation (10). In the integration by parts of Equation 
(10) the boundary terms have been neglected. Next, the trans- 
formation r, = r cos <p and - r sin <p gives rise to two terms, 
one of which is Equation (11), and the other has been neglected. 
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2.5 — Continued. 

If the correlation distance T is small compared to the diameter 

of the surface element 6 S, then the factor X 2 ( v z * ~ V z ) 

_ ^ x* ( v ) will be small except near r = 0. This fact 

makes the neglected term small. 

e. The scattering formulation presented here considers only the 

amplitudes of the incident and reflected fields; thus, no polari- 
zation information about the scattered field is available. The 
only polarization information comes from the Fresnel reflection 
coefficients. A vector theory for scattering of a plane wave 

from a plpna surface with random irregularities has been dis- 

[4 5] 

cussed in numerous articles. 1 

[3] 

£. Shadowing has been taken into account as suggested by Beckman; 

this has been discussed in Section 2. 3. 

g. Multiple scattering has not been considered. 

h. While the basic theory of Section 2. 1 assumes a plane wave incident 
on a plane reflector, we have modeled a point source incident on a 
sphere by breaking up the surface of the sphere into elements 
which are approximated by plane sections. The field incident on 
each section may then be approximated by a plane wave. The 
resultant reflections are added up to yield the total field at the 
receiver. Some authors 11, 6J have integrated over a plane earth 
and modified the result by the divergence coefficient; we call this 
the flat earth approximation. This procedures does not satisfy our 
intuition, and there are other difficulties associated with it. First, 
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2.5 — Continued. 


the region of integration is completely unspecified. Second, the 
far field approximation is violated. We have attempted to calcu- 
late the received field using this approximation; the differences 
in the results are dramatic in most cases, as indicated by 
Figure 3-1. 

i. The three-fold integral Equation (11) has been computed numerically. 
In Beckman^ approximations are obtained for small and large 
values of g, where g is the roughness criterion. These approxima- 
tions are not valid for moderate values of g, and since all values 

of g are liable to arise, we have not used these approximations. 

j . The surface height function is assumed to be normally distributed 
Equation (1). 

2.6 Pulse Delay and Broadening. 


To evaluate the effect of multipath on a communication system, it is necessary to know 
the time delay between the multipath and the direct path and to know to what extent a 
pulse is broadened in time by reflection from the earth's surface. The time delay 
between the direct and multipath pulses is easily calculated from the geometry of a given 
problem. In Figure 2-4 a typical geometry is shown. The time delay is given by 

R + S - Q 

T = 

C 


where C is the velocity of light. 
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2.6 — Continued. 

To calculate pulse broadening it is necessary to divide the scattering surface into 
annular regions, each of which has a nearly constant time delay. To calculate the 
received field, the surface integral in Equation (13) is evaluated over this annular 
region. 
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3. RESULTS OF CALCULATIONS. 


In this section some results of calculations using the theory of Section 2 are discussed. 
The only effects included in these calculations are rough surface scattering and space 
loss. Pulse distortion, Faraday rotation, etc. , are not included in the present report 
results. All are for VHF frequencies and for a terrain corresponding to marshy land. 

3. 1 Discussion of Results. 

f fll 

Figure 3-1 shows a comparison between the RCA model 1 J and our curved earth and 
flat earth models. Our calculations are performed for a roughness, o, of 1. 0 meters 
and a correlation distance, T, of 10 meters. In the RCA model the ration o/T is 
assumed small and does not enter explicitely into the calculation. The altitude of one 
satellite is 275 nautical miles and the other is at synchronous altitude. Our curved 
earth model shows a dip at around 45 degrees which neither of the other two calculations 
show. This dip is caused by the fact that the diffuso scattering is decreasing while the 
specular term is increasing as the incident angle increases. Our curves approach 
zero at 90 degrees due to the divergence term, D, and the shadowing term, 8(0^. 

The RCA curve approaches zero due to the divergence term and due to the fact that 
it does not contain the specular term. 

Figure 3-2 shows the effect of varying the altitude of one satellite while the other stays 
at 275 nautical miles. Roughness is 1.0 meter and correlation length is 10 meters for 
all calculations. 

Figure 3-3 shows the effects of varying the roughness and correlation length. One 
satellite is at 275 nautical miles and the other is a synchronous satellite. For one 
value of roughness, 1 meter, the minimum in the multipath signal moves to smaller 
incidental angles as the correlation length is increased from 10 motors to 30 motors. 
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INCIDENT ANGLE AT SPECULAR POINT (DEGREES) 

Figure 3-1. Multipath Signals Using Three Different Calculation Schemes 


CO 

I 

to 


ESL-PR53 






to 

I 

to 


Figure 3-2. Multipath Return for Satellites at Various Altitudes (in Meters) 
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3. 1 — Continued. 

It is also seen tha t for one value of the correlation length, 10 meters, the minimum in 
the returned signal becomes less pronounced as the roughness goes from 1 meter to 
4 meters. It is apparent that there is considerable dependence upon a and T in the 
multipath returned signal. 

Figure 3-4 shows the returned signal as a function of time for one satellite at 275 nautical 
miles and the other satellite at synchronous altitude. A pulse duration time of 1 micro- 
second was used and the returned signal was calculated in 2 microsecond steps out to 
12 microseconds. The curves have a maximum in the zero to 2 microsecond region 
because of the specular signal. Beyond 2 microseconds the signal stays very constant. 


3.2 Conclusions. 

From Figure 3-1 it is apparent that there are considerable differences in the resulting 
multipath returns between the flat earth approximation and the more appropriate 
spherical earth case. This is largely due to the poorly defined scattering area for the 
fla t, earth approximation. 

The calculation presented in Figure 3-2 shows that there is some dependence of the 
multipath return on the altitude of the two satellites. This would be expected because 
of the change in the size of the scattering area and the change in the space loss as the 
altitudes vary. 


Figure 3-3 is quite important because it shows for the first time the exact dependence of 
the multipath signal on the parameters o and T. It can be concluded that there is a defi- 
nite dependence on a and T and that it should not be ignored in evaluating the effect of 
multipath on a satellite to satellite communications system. 
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3.2 — Continued. 


The pulse versus time calculations show that it is necessary to perform this type of 
calculation to determine at what satellite position the most multipath return will be 
received. It can not be concluded that the most serious interference case for a pulsed 
system occurs when the satellites are overhead. For a CW tr ansm itter the worst 
case will be dependent upon data rate and pulse duration time, but will occur when the 
two satellites are directly overhead. 

A final conclusion which was reached during the derivation of the theory is that previous 
calculations based on approximations to Equation (11) are of doubtful value. The reason 
for this is that, the approximations are valid only for the case of very diffuse scattering. 
As either the roughness gets smaller or the separation between satellites gets greater, 
these approximations are very suspect. 
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